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The convective instability in a plane liquid layer with time-dependent temperature 
profile is investigated by means of a general method suitable for linear stability analysis 
of an unsteady basic flow. The method is based on a non-normal approach, and predicts 
the onset of instability, critical wavenumber and time. The method is applied to transient 
Rayleigh-Benard-Marangoni convection due to cooling by evaporation. Numerical results 
as well as theoretical scalings for the critical parameters as function of the Biot number 
are presented for the limiting cases of purely buoyancy-driven and purely surface-tension- 
driven convection. Critical parameters from calculations are in good agreement with those 
from experiments on drying polymer solutions, where the surface cooling is induced by 
solvent evaporation. 



1. Introduction 

Thermally driven flows near liquid interfaces continue to be an active area of research 
since the first experimental studies by H. Benard over 100 years ago. They are character- 
ized by a multitude of interacting physical mechanisms and display a large variety of reg- 
ular and complex flow patterns. Experimental and theoretical studies of such flows have 
stimulated and accompanied the development of the theory of pattern formation and non- 
linear phenomena in general. Recent develo pments are summariz e d in a number of mono- 



graph s and review papers (see fo r example IColinet et al\ (|2001l ); iNepomnvashchv et al 
l|2006t ); lBodenschatz et all (|2000h ). 



The driving forces in these thermally driven flows are surface-tension gradients (Marangoni 
forces) and buoyancy forces, which originate from the temperature dependency of surface 
tension and density, respectively. Both mechanisms are classically studied in convection 
taking place within a finite layer subjected to a steady temperature difference: for the 
buoyancy - drive n case, this is the celebrated Rayleigh-Benard convection formulated by 
Rayleigh (1916), for the surface-tension driven case this is the Benard-Marangoni con- 
vection studied by iPearso 3 <|l958h . In both CcLSGS, cL certain critical temperature differ- 
ence must be applied in order to sustain convective motion. The theoretical predictions 
for such critical temp e rature differences a re in excellent agreement with experiments 
( Chandrasekhar ( 1961 ): Schatz et al. ( 19951 )) . Thermal convection can also appear when 
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the external conditions are time-dependent, e.g. by a modulation or abrupt change of 
cooling or heating. In this case, the basic conductive temperature distribution is time 
dependent, and the stability problem becomes a non-autonomous one, i.e. the solution 
cannot in general be sought in the form of exponentials exp (crt). For time-period i c mod - 
ulation of the basic s t ate on e can resort to Floquet theory (jRosenblat fe Tanakal (|197lf ); 
Bhadauria fe Bathial (J2002)) but, in the general case, the formulation and prediction of 



critical conditions for the onset of convection becomes much less clear-cut than for a 
steady base state. The following two basic approaches are common for a general, time- 
dependent basic state: 

• reduction to an autonomous problem by the frozen-time assumption, whereby the 
basic state is supposed to evolve much more slowly than the perturbations, and 

• solution of the full non-autonomous linear perturbation problem for some initial 
conditions which are supposed to be representative. This is called amplification theory 
(denoted by AT in the following). 

Both approac hes go at least back to the 1960s . The first one w as used b y Lick (196^) and 
Currid (1967) and the second by iFoster (1965). Later work by iHomsvl ( 1973 ) introduced 



energy stability ideas, but the resulting bounds are not necessarily useful for predicting 
instability thresholds and the mathematics is considerably more involved. In each of these 
works the focus was on buoyancy-driven convection. 

The frozen-time approach was used for the Marangoni convection bv lKang fc Choil (|1997l ). 
These authors studied the dynamics of a fluid layer subjected to a sudden change in sur- 
face temperature. However the frozen-time method was only applied for late times in the 
evolution. For early times they formulated the so-called propagation theory, whereby the 
problem is again reduced to an autonomous one by using a certain similarity solution for 
base state and perturbations. This approximation is actually appropriate for an infinite 
layer. The frozen-time model relies on the validity of the assumption of fast instability de- 
velopment (relative to the base state). This could be inadequate at instability thresholds, 
at which the instability can develop on the same time scales as the base state. The other 
method namely the amplification model, has a sounder mathematical foundation, but 
comes with the necessity of identifying representative initial conditions and amplification 
levels for the perturbations. 

A specific analysis suitable for transient problems is developed in this paper, which is 
an extension of the previous amplification model. The underlying concept is called the 
non-normal approach. In such a case, the choice of initial perturbations is not some- 
what arbitrary as in AT but it is guided by an optimization technique. This procedure 
consists in identifying, for a given time t and a given wavenumber k, the initial per- 
turbation with the strongest amplification (called the optimal perturbation). Hence this 
method may exhibit initial conditions which are not a priori obvious configurations. As 
a consequence it gives the maximal value for the perturbation norm that was previously 
defined for the problem under study. This cannot be achieved through an AT simula- 
tion. Concerning the determination of initial conditions, we note that the non-normal 
approach has also been used for insta bility problems of steady base flow, e.g. in parallel 
shear flows, such a s boun dary layers (jSchmid fe Henningson (2001)) or plane Poiscuillc 
flow ( Reddy et all ( 19981 )). For such problems, there was a long-standing gap between 
the standard two-dimensional Tollmien-Schlichting modes obtained using the classical 
normal mode analysis - though these modes were only observed in some very controlled 
experiments - and the general experimental observation of streak generation. It is by 
resorting to non-normal mode analysis that the coherent structures observed in bound- 
ary layer transition, i.e. streamwise rolls and streaks, could be deduced from an optimal 



3 



approach and their generatio n linked to the lift-up mechanism proposed by Landahl 
(jSchmid fc Henningsonl <|200lh >. 

The non-normal approach can be extended to transient cases. This is precisely what 
we propose here. Such a method provides the strongest amplification, i.e. the optimal 
growth at a given duration T after initial conditions have started their evolution. If this 
amplification is sufficiently large, the optimal perturbations may result in the modifi- 
cation of the basic flow at time T and we can speak of an unstable regime in a way 
to be defined by a norm. To the best of our knowledge, this approach is the only one 
capable to provide clear-cut answers on instability problems for truly unsteady basic 
flows. However, as mentioned by one referee, this definition of the stability condition is 
blurred owing to a certain arbitrariness in choosing the norm and defining the critical 
amplification gain. This is inherent to any unsteady linear stability problem and does 
not depend on the particular method. It is then more suitable to view the results as the 
estimation of a transition region between a domain exhibiting strong convection and a 
domain where initial perturbations are damped or have no time to significantly develop 
during the transient regime. The non-normal approach allows one to characterize this 
transition domain properly. The results presented here for different norms and different 
critical amplification values show that this transition region is thin, so that the notion of 
stability threshold is still valid for this transient problem. On the contrary, the frozen- 
time assumption may fail, for instance if the base state evolves on the same time scale or 
faster then the unstable modes characterizing its frozen-time stability. In that case, the 
computed growth rates might not correspond to any true amplifications, at the system 
time scale. This might affect the determination of critical conditions 0- 

In the present work, we provide the results obtained by means of the non-normal 
approach as well as the frozen-time approach. In the specific flow case presented here, 
the quasi-static method is shown to provide similar results except for some particular 
quantities. 

The proposed analysis is general and can be easily extended to many oth e r unst eady 
problems, e.g., chemically driven hydrodynamic instabilities (jEckert et al. (|2004 )). In 



the present paper, we determine the onset of convection in a drying experiment, which 
leads to an unsteady Benard-Marangoni problem. More specifically we study the sudden 
cooling due to evaporation of a liquid layer, where the decrease of surface temperature 
is induced by the vaporization latent heat. This pro blem has been th e subject of many 



experi mental or theoretic al studies (see for examp le iBerg et all |Vidal fe Acrivosl 

(|l968l ), and more recently [Mancini fc Mazal (|2004f) or lMoussv et all (12004)). The specific 



motiv ation for our work is the drying process studied in experiments bv lToussaint et 



(2008) aerformed on a polymer solution (Polyisobutylene/Tolucne). The solution initially 



at the ambient temperature is poured in a dish located in an extractive hood. When evap- 
oration of toluene begins, convective patterns are observed at the very beginning of the 
experiment (quasi-instantaneous or less than 100 s after pouring the solution). They dis- 
appear well before the end of the drying. The very large Lewis number Le = n/D mo i 
(k and D mo i denote the thermal and mass diffusivity) of the polymer solution (about 
10 3 ) is an indication that the thermal diffusion is faster and that convective patterns 
observed in the first minutes should be mainly driven by thermal effects. Two experi- 

t In the frozen-time assumption, the critical conditions are determined using the quasi-static 
growth rates cr(t). However the pertinent character of this method stongly depends on the 
rapid variation of these "quasi-static growth rates" with the basic state changes. For instance, 
the amplification between times ti and ti would be, at zeroth order WKB theory, equal to 
f^ 2 a(t)dt. So if the instantaneous growth rate strongly varies during this interval, its positive 
value at a given moment may not be interpreted in the right way. 
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mental observations detailed in iToussaint et al. (|2008l ) support this thesis. First, a few 
experiments were conducted with deuterated solvent, whose density is higher than the 
polymer density. In that case, the density of the solution decreases when the polymer 
concentration increases, leading to a stable situation if the solutal Rayleigh-Benard prob- 
lem is considered. Since no differences were found with the experiments conducted with 
the standard solvent, we can exclude solutal buoyancy as a dominant mechanism. Sec- 
ond, free surface temperature fields measured by infrared camera showed that the end 
of free convection was related to the duration of the transient thermal regime. In these 
free convection experiments it can be inferred from the work on steady convection by 
Pearson ( 1958f ). that the Marangoni effect is dominant for thicknesses typically less than 



lem and the buoyancy dominant for higher th icknesses. For a more accurate description 
of the experiments see Toussaint et al. ( 20081 ). 



The paper is organized as follows. In section [2l we present the basic assumptions of 
the model and the governing equations. Thereafter the unsteady basic state is described 
and a specific stability analysis is introduced in section [3l In particular the non-normal 
method is explained and the choice of norms is discussed. Section 0] contains the main 
numerical findings for the two limiting situations: the pure Rayleigh-Benard and the pure 
Marangoni case. Critical conditions for the optimal modes are presented. Part of these 
numerical results are also obtained by a scaling analysis. A comparison of this method 
with the frozen time approach is also performed. Finally, in section [5j the comparison 
with experimental results is discussed. 



2. Mathematical Model 

2.1. Basic Assumptions of the Model 

The mathematical model of Rayleigh-Benard-Marangoni convection used throughout this 
paper is based on a one- layer model in which three assumptions are made : (1) the upper 
surface remains planar, (2) the layer thickness d remains constant, (3) the heat and mass 
fluxes across the upper surface are given by transfer coefficients. Moreover our analysis 
is restricted to fluids characterized by a Prandtl number Pr = v/n ^ 1 with v the 
kinematic viscosity. This turns out to be the case of most liquids (including water and 
organic solvents). In that instance, the thermal diffusion time scale is always larger than 
the viscous diffusion time scale. 

The first hypothesis can be tested as f ollows. Two modes o f instability are known to oc- 
cur in t he Benard-Marangoni pr oblem ( Scriven fe Sternlrngl ( 1964 ) ; Reichenbach fc Lindei 



(jl98lh : iGoussis fc Kelhl (jl990h V One is generated from the interaction of the velocity 



perturbations with the basic temperature, while the other mode, characterized by a wave- 
length long compared with the layer thickness, is due to the coupling of the Marangoni 
effect with the deflection of the free surface. Mathematically, surface deformation can be 
neglected on scales of order d when the Laplace pressure asso ciated with a curvature 1/d 
is large compared with the dynamic pressure. This condition ( David ( 198?t )) corresponds 



with the smallness of the crispation number Cr = (pvn) / (ad) where p and a respectively 
denote fluid density and surface tension. Moreover, the Galileo number Ga = (gd 3 ) / '(vk) 
characterizes the relative importance of gravity (g gravity constant) and diffusion. A 
large value of the Galileo number indicates that gravity stabilizes the long-wave mode. 
The free surface deformation can be neglected if Cr <C 1 and Ga ^S> 1. Such conditions 
are shown be satisfied for the experiments considered here. 

The second hypothesis can be adopted if Pe <C 1, where the Peclet number Pe = (dv ev )/ K 
is defined as the ratio between v ev the interface velocity due to evaporation which is equal 
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to minus the time derivative of d(t), and the thermal velocity scale n/d. Indeed when 
Pe <C 1, the surface displacement v e v5diff remains negligible compared to the total 
thickness d during the problem characteristic time i.e. the diffusion time Sdiff = d 2 j k. 
In the experiment (see section[5]), the Peclet number is smaller than 0.1. 
Finally let us discuss the third assumption. The boundary condition at the free surface 
results from the coupling between the system and its surroundings. In evaporation ex- 
periments, the solvent flux and thus the temperature gradient in the fluid depends on the 
heat and mass transfer with the ambient air. Several authors have de veloped numerical or 



theore ti cal studies tak i ng int o account this coupling, see fo r examp le Mcrkt & Bestehorn 
( 20031) ; IColinet et cd\ (|2003l) ; lOzen fc Naravanar] (|2004h ; iMoussv et all (|2004 . In this 



paper we adopt a simple description by global heat and mass transfer coefficients, as our 
main interest is directed on the transient character of the problem under study and not 
on the detailed description of the transfer per se. 

2.2. Governing equations 

We formulate the basic equations in a Cartesian coordinate system, where the bottom 
of the layer coincides with the plane z = and the upper free surface with z = d. The 
fluid is characterized by a density p, a kinematic viscosity v 1 a thermal diffusivity K and 
a thermal expansion coefficient a. The Boussinesq approximation is assumed to govern 
the velocity field v — v x e x + v y e y + v z e z and temperature field T 

^ + (v-V)v = -^ + is\7 2 v + ga(T-T oc )e z , V ■ v = 0, (2.1) 
at p 

dT 

— + (v-V)T = K V 2 T, (2.2) 
at 

where g denotes the acceleration due to gravity and the temperature of the ambient 
air. In this approximation, the density p is taken to be the density at T = T^. At the 
bottom, the velocity satisfies the no-slip condition and the wall is assumed adiabatic 
since, in the experiment, the bottom of the dish is thermally insulated by an air gap. 

v = 0, d z T = at z = 0. (2.3) 

The upper boundary conditions are more involved. Assuming a planar surface, the local 
evaporation mass flux reads 

Q m (T s (x,d,t)) = p(v z (x,d,t) - j t (d)) (2.4) 

Moreover the global mass balance reads (no fluid is introduced to counterbalance the 
evaporation mass loss) 

Qrn = -P^W = P V ev (2-5) 

dt 

where Q m is the mean evaporation flux over the free surface. From relations 12.41 and 12.51 

we get : 

v z = pr v ev (2.6) 

If we assume that the flux variations Q m (T s (x 1 d,t)) — Q m are much smaller than the 
flux Q m itself, then \v z \ v ev = Pe n/d, the Peclet number being defined using the 
evaporation velocity (see section 2.1). Since Pe <§; 1 and njd is the velocity scale, the 
kinematic boundary condition so reduces to 

v z = at z = d. (2.7) 
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In addition, the balance of tangential forces at the upper surface requires that the velocity 
field satisfies 

pvd z Vj = djd{T s ), (j = x,y) &tz = d, (2.8) 

where T s denotes the temperature inside the fluid layer at z = d and cr(T) denotes the 
surface tension which is assumed to be a linearly decreasing function of temperature, 

*(T) = *(T 00 )- 1 (T-T 00 ), 7 = -^>0. (2.9) 

Finally, the conservation of energy flux should be imposed at the interface. It reads : 

dT 

-\-^ + h th {T QO -T s )=LQ m {T s ) (2.10) 

The first l.h.s. term represents the heat conduction in the liquid, A denoting the thermal 
conductivity related to the thermal diffusivity via A = npC with C the liquid specific 
heat. The second l.h.s. term expresses the heat flux density in the gas using a simple 
phcnomenological model based on the heat transfer coefficient h t h and the difference 
between the air temperature far from the liquid and the surface temperature T s . 
Finally the cooling effect due to solvent vaporization is expressed in the r.h.s. term, 
where L stands for the latent heat of vaporization and Q m for the solvent mass flux per 
unit area. This latter quantity depends on the surface temperature T s . In the framework 
of the one-layer model, one imposes 

Q m {T s ) = h m (c s (T s ) - Coo) (2.11) 

where c s (resp. Coo) denotes the solvent concentration in the gas phase near the surface 
(resp. far from the surface) and h m is the phenomenological mass transfer coefficient in 
the gas. Assuming local thermodynamic equilibrium, c s directly depends on the surface 
temperature through the saturated vapour pressure. The variable Q m can be linearized 
around which leads to the final expression: 

dT dO 
-X—+H th (T 0o -T s ) = LQ m {T oo ), H th = h th + L-^{T 00 ). (2.12) 

Initially the liquid layer is isothermal i.e. T(z,t = 0) = T^. At large times, the system 
reaches a steady state in which the temperature in the layer is again uniform but with 
a temperature difference — T s — AT st > with the gas located far from the surface. 
This difference is imposed by the condition (|2.12|) . where 

AT st = (2.13) 

To put the above equations in a non-dimensional form, scales for temperature, length, 
velocity and time are needed. The temperature difference AT st provides the temperature 
scale and the layer thickness d the relevant length scale. Two velocity scales can be intro- 
duced in this problem namely the viscous velocity scale V = v j d and the thermal velocity 
scale V = n/d. Here only the thermal scale V = n/d is used. Finally the appropriate 
time scale is the thermal diffusion time d/V = d 2 / k. The non-dimensionalization leads 
to equations for v and 6{z,t) = (T(z,t) - T 00 )/AT st : 

^-{d t v + (v ■ V)v) = -Vp + V 2 v + Ra9e z , (2.14) 
Pr 

V-v = 0, (2.15) 

d t + (v ■ V)0 = V 2 6», (2.16) 
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d z v x + Mad x 9 = d z v y + Mad y 6 = at z = 1. (2.17) 
d z 9 + Bi9 + Bi = at z = 1, <9 2 <9 = at z = 0, (2.18) 
v z = 0, at z = 1 u x = w y = ^ = at z = 0. (2-19) 

in which the Rayleigh, Marangoni, Prandtl and Biot numbers 

Ra=^^i, Ma,*, Pr= V -, = ^ (2.20) 

appear. 

In the following, we discriminate between two opposite cases: Bi < 1 or 1 < Bi. This 
corresponds respectively to an effective conductance H th in the gas much smaller or much 
larger than the heat conductance X/d in the fluid. 



3. Basic state and optimal linear perturbations 

We study the stability of a purely conductive unsteady basic state 9q(z, t) which is initially 
uniform i.e. 6o(z,t = 0) = 0. This state evolves since the upper free surface is cooled 
down by the latent heat released through evaporation : The velocity field remains always 
zero but the unsteady field 9q(z, t) satisfies a pure heat equation 

pin 3^0 

— - = — -£ with d z 9 = at z=0, d z 8 + Bi6 + Bi = at z=l (3.1) 
at az 1 

The basic state only depends on the Biot number Bi. Note that the term Bi9 character- 
izes the heat transfer in the gas phase and the term Bi represents the cooling effect due 
to evaporation. The unsteady temperature field 9o(z,t) is shown in figure [1] for different 
Biot numbers and times. A cooled layer develops from the upper surface whose thickness 
^o(i) ~ min(\/t, 1) is similar at a given time for different Biot numbers and increases 
until it fills the whole layer. Conversely, if A9o(t) denotes the characteristic temperature 
difference within the fluid, the maximum reached over the whole time evolution by this 
quantity increases with Biot number Bi (figure [TJi) . Actually, two different regimes are 
observed according to the value of Bi\^. For small Biot numbers, typically Bi < 1, the 
cooled layer reaches the bottom while the jump A6*n(i) is less than one. Afterwards, 
A6o(t) decreases. Such time evolution can be summarized by the following scalings 

|0 O | ~ A0 O ~ BiVi, S (t)^Vt for < y/i < 1, (3.2) 

|0 o |~-Bi*, A0 Q (t)~Bi, 5 (t)~l for l<t<Br\ (3.3) 

For Bi -1 < t, the temperature field #o relaxes towards the steady uniform temperature 
9q(z, t) = —1 and all the energy needed by evaporation is carried on by convection in the 
gas phase. 

For large Biot numbers, typically Bi > 10, the temperature jump A0o(t) reaches a max- 
imum before the cooled layer reaches the bottom (figure [TJi). In that case, the maximal 
jump is equal to one and \9q(z = l,t)\ ~ 1 at that time. The time evolution can be 
summarized by the following scalings : 

|0 O | ~ A0 O ~ BiVi, S (t)~Vi for < t < Br 2 , (3.4) 

A6» ~ 1, S - Vi for time Bi~ 2 < t < 1. (3.5) 

For 1 < t, the temperature decreases in the whole layer thickness to reach the steady 
state regime 9 (z, t) = —1. 



t Details are contained in Appendix A. 
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In order to study the stability of the unsteady conductive state we split the fields into a 
basic flow and three-dimensional perturbations 



v p (x,y,z,t), 



h(z, t) + 6 p (x, y, z, t), p = p (z, t) + p p (x, y, z, t) 



(3.6) 



and linearize in the perturbations to get a set of linear equations. Since this problem 
has no preferential direction in the (x, y) plane, the perturbation Fourier modes in such 
directions decouple in the linear regime. Without loss of generality, we thus consider a 
nondimensional wavenumber k in the x direction and no dependence in the y direction 
reducing the flow to be two-dimensional 

(v p ,6 p ,pp) = (u(z,t),0,w(z,t),9(z,t),p(z,t))exp(ikx) 

These infinitesimal perturbations are governed by the linear system 

d 2 



19. ., . 
Pt8i U + tkp - 



dz'- 



- k< 



J_d_ 
Prdt 



dp 

dz 



k 



w = 0, d z u 



d?_ 

dz 2 

dt 

Ma ik6 = 0. 



Ra9 = 0, 

dz 
dj 



, A dw 
iku + — — 

dz 



dz- 



-k 2 



w 



Bi9 = 0, at z - 
at z 



0, ^ 

dz 





(3.7) 


o, 


(3.8) 


o, 


(3.9) 


o, 


(3.10) 


1, 


(3.11) 


0. 


(3.12) 



To quantify the amplification gain at time t\, it is customary to define a norm which is 
generally based on the kinetic energy of perturbations. This norm can be orthogonally 
decomposed in a Fourier basis so that the individual contributions of each wavenumber k 
can be studied independently. In the present case, the temperature field is playing a 
major role as well. We thus define two different norms corresponding to two different 
situations. The first one is based on the kinetic energy of perturbations 



Ev(h) = / (u(z, ii)u + (z, ti) + w(z, ti)w + (z, ti))dz 



(3.13) 



where superscript + denotes complex conjugation. The integration is performed over 
the entire layer width and perturbations are obtained after integrating the above linear 
system over the time period [0, <i]. The second one 



E T (h)= / 9(z,t 1 )9+(z,t 1 )dz 



(3.14) 



is based on the temperature field and not the velocity field. 

For finite Prandtl number two extreme cases are considered, with the initial perturbation 
concerning either the velocity field, or the temperature field. In the following we will 
consider the amplification factors Ev(t\)/Ev(0) or Exitx)/ Et(0) to characterize the 
stability. Then, when the initial perturbations only concern the velocity field, we only 
take into account the amplification Ey(ti)/Ey(0) since £t(0) = 0. Conversely, when the 
initial perturbations only concern the temperature field we use the amplification factor 
£t(£i) j 'Et(0). For the infinite Prandtl number case, the velocity perturbations are not 
dynamical quantities since the time derivative of the velocity drops out from equations 
(|3.8p - (|3.9p . In other words, the velocity perturbations are slaved to the temperature 
perturbations. In this case, we use only £ , t(^i)/^t(0) and hence the norm Ej 1 . 
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For a given initial disturbance, one evaluates the amplification gain at time ti by com- 
puting E(t\) / E(Q) where E — Ev or E — Et- It is the purpose of a non-normal 
analysis to compute the quantity G(t±; k; Ma, Ra, Bi, Pr) = Max[E(ti)/E(0)] which 
is the upper bound for the energy amplification that a disturbance of wavenumber k 
can reach at time t±. This approach solves a sort of the finite time stability problem : 
it investigates the transient evolution and defines for any time t-\, an optima l pertur- 



bation mode which actually reaches the upper bound (jSchmid &: Henningsonl ( 200ll )). 
This mode (i.e. the optimal initial z— profil e) is found numerically by solving an op - 
timization problem dFarrell fc Ioannou (|1996h : lAndersson et al. ni (2000): 
Schmid fc Henningsonl ( |200ll )). The optimum of E(t\) is determined taking into account 
several constraints: (i) the disturbance energy at time t = is equal to unity; (ii) the 
disturbance satisfies the linear governing equation as well as the boundary conditions 
during the complete time interval [0,ti]. This problem is best solved with the help of 
a Lagrangian formalism and Lagrangian multipliers which are introduced to precisely 
enforce the above constraints. In the present case, these multipliers are adjoint fields 
(u(z, t), w(z, t), 9{z, t),p(z, t)). Following a standard derivation, these quantities satisfy a 
set of adjoint equations. It is 



if 



k 1 



1 d _ dp 
Pr dr dz 



dz 2 



-90o 
dz 



dz 2 
0, iku + 



u = 0, 
dw 



OT 



8?_ 

dz 2 



m 



dz 

— Raw = 



u = w = 0, — = at z = 0, 
oz 

w = 0, d z il = 0, dJ + Bi6 + Mad z w = at z = 1, 



(3.15) 

(3.16) 

(3.17) 

(3.18) 
(3.19) 



in which r = —t. These adjoint equations have to be solved backwards in time. Let us 
denote by the symbol q the vector field (u, w, 0,p). One obtains the optimal perturbation 
for time t\ by an iterative scheme which propagates a given initial condition forward in 
time using the direct problem (here denoted by Fj(q) = 0, j = 1 . . . 4), the result of which 
serves as an "initial" condition for the backward propagation by the adjoint equations 
(here denoted by Fj(q(t)) = 0, j = 1 . . .4). More specifically Q, a relation between the 
adjoint q(z, t\) and q(z, t\) is imposed. After one forward-backward integration, quantity 
q(z, 0) is obtained and a relation between q(z, 0) and q(z, 0) is also imposed. An updated 
initial condition for the next iterative step is then available. This process should be 
self-consistent : one uses an iteration procedure which is schematically illustrated by a 
diagram 



q(z,0) 



F,(q)=0 



Fj(q)=0 



q(z,h) 

i 



(3.20) 



q(z,0) < — q(z,ti) 

Convergence is reached when the initial condition for the forward problem does not 
change appreciably - up to a normalization constant - by an appropriately chosen cri- 
terion from one iterative step to the next. The converged mode is precisely the initial 
optimal perturbation for time t\. The maximum energy amplification is computed by 



f Details are contained in Appendix B. 
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propagating the converged initial condition forward in time and by forming the ratio of 
the disturbance energy at the end of the time interval to the energy at the beginning. 
The direct and adjoint equations have been discretized using a pseudospectral scheme 
based on Chebyshev polynomials and a streamfunction-based formulation to account for 
incompressibilitjU 

4. Numerical Results 

4.1. Quantities provided by the non-modal analysis 

For the unsteady basic state 9o{z,t), a given set (Ma, Ra, Bi, Pr), and a given type 
of initial perturbation (temperature or velocity), the non- modal analysis determines at 
a given time t±, the maximum energy amplification G(t\; k; Ma, Ra, Bi, Pr) over all 
possible perturbations of wavenumber k (see figured]). In a way, the value InG/ii is 
analogous to a growth rate for classical stability analysis. More generally, it appears 
possible to extend the usual concepts of classical stability analysis to unsteady flows. 
For instance, G(ti, k; Ma, Ra, Bi, Pr) can be maximized over wavenumber k and time t\ 
providing the global maximum amplification G max (Ma, Ra, Bi, Pr) (see figure [3]) ■ This 
value is reached at time t = t opt (Ma, Ra, Bi, Pr), for an optimal wavenumber denoted 
by k op t(Ma, Ra, Bi, Pr) and for a specific perturbation structure in z. These latter two 
quantities play the role of the most amplified wavenumber and of the most amplified 
mode for the standard stability analysis. One also obtains a "stability" diagram, by de- 
termining the region of the space (Bi, Pr, Ma, Ra) where the amplification gets above 
a threshold Gthres- It is a way of separating the region where amplification or attenu- 
ation occurs. The value of Gthres-, for instance Gthres = 1, is somewhat arbitrary : as 
already said in the introduction, an unsteady problem is indeed characterized more by a 
transition domain than a well defined threshold. It is demonstrated below that choosing 
Gthres — 1 or Gthres = 100 does not result in major differences, so that the transition 
region is thin compared with the absolute values of the Marangoni and Rayleigh num- 
bers. For fixed Ra and Pr, one can determine the curve Ma c (Bi, Ra, Pr) such that if 
Ma < Ma c (resp. Ma > Ma c ) then G max < Gthres (resp. G max > Gthres)- Similarly 
one may define for fixed Ma and Pr, the curve Ra c (Bi, Ma, Pr). These curves play a 
role very much similar to marginal stability curves. Each point of the critical curve is 
associated to a critical wavenumber k c = k op t and critical optimal time t c = t op t. Note 
that until this point, most of this procedure can be extended to other unsteady problems 
in a straightforward way. A comparison between these critical curves and the experimen- 
tal diagram that separates the domains where convection is observed or not observed, is 
made in section [5] 

4.2. Infinite Prandtl number: the pure Marangoni case Ma ^ and Ra = 0. 

For infinite Prandtl number, the velocity is slaved to the temperature field. As a conse- 
quence, only perturbations in temperature field are pertinent. In the plane (Bi, Ma), the 
critical curve Ma c (Bi), critical wavenumber k c (Bi) and critical optimal time t c (Bi) are 
presented for the pure Marangoni case and two thresholds Gthres = 1 arid Gthres = 100 
(figured]). The critical Marangoni number Ma c (Bi) slightly depends on the value of the 
threshold and seems to be consistent, within the numerical uncertainties, to the following 
laws (here Gthres = 1) 

Ma c (Bi) ~83/Bi for Bi < 1 and Ma c (Bi) ~15Bi for 1< Bi (4.1) 



£ Details are contained in Appendix C. 
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Moreover the wavenumber k c {Bi) (figured) is an increasing function of the Biot number 
while optimal time t c {Bi) is a decreasing function of the same parameter. 
Using heuristic arguments 0, the following scaling laws can be deduced : 
For Bi -C 1 the critical conditions for instability can be expressed as 

Ma c ~\/Bi with V~Bi < k c < 1 and \<t c <l/Bi. (4.2) 

For 1 <C Bi, the critical conditions become 

Ma c ~Bi with 1 < k c < Bi and t c ~ k~ 2 (4.3) 

The spatial structure in z of the optimal perturbation at k — k c and Ma — Ma c is 
shown on figure [5] for two Biot numbers and two different times : time t = and time 
t = t c when the perturbation reaches its maximum amplification. The spatial structure 
of this optimal perturbation is shown to change slightly during the time evolution. In 
this respect, this optimal mode does not differ much from the classical most amplified 
mode of steady problems. 

4.3. Infinite Prandtl number : the pure Rayleigh case Ra ^ and Ma = 0. 

In the plane (Bi,Ra), the critical curves Ra c {Bi), k c {Bi) and t c (Bi) are presented for 
two thresholds Gthres — 1 and Gthres — 100 (figure [6]). The critical Rayleigh Ra c (Bi) 
slightly depends on the value of the threshold and seem to be consistent, within numerical 
uncertainties, to the following laws (here Gthres = 1) : 

Ra c (Bi) ~ 600/ Bi for Bi < 1 and Ra c (Bi) ~ 960 for 1 < Bi (4.4) 

Moreover the wavenumber k c (Bi) (resp. time t c (Bi)) is an increasing (resp. decreasing) 
function of the Biot number for small Biot and reaches a plateau for larger Biot number. 
One can deduce using heuristic arguments 0, the following scaling laws: 

Ra c ~ l/Bi with \J~Bi < k c < 1 and 1 < t c < l/Bi, for Bi -C 1 (4.5) 
Ra c — 1 with k c ~ 1 and t c ~ 1, for 1 < Bi. (4.6) 

4.4. Comparison with classical steady results and with transient frozen-time method 

It is wor t h co mparing t he results presented h ere with the well-known results obtained by 
iPearsonl (|l958l ) § (resp. I Sparrow et all (|l963l )) in the framework of the steady Marangoni 



(resp. Rayleigh) problem. This comparison is pertinent since the boundary conditions at 
the top and the bottom of the layer (equations 13. 11113. 12"]) are similar in the present 
work and in these classical analyses. However, the Marangoni and Rayleigh numbers 
defined by these authors are based on the steady temperature difference ATq between 
the top and the bottom of the layer. This steady temperature difference is missing in 
the transient problem under study. It is then not possible to make a direct comparison 
between the thresholds values obtained in our paper and the previous ones from Pearson's 
or Sparrow's publications, and a preliminary transformation is needed. Indeed, at each 
time t, one might define the equivalent temperature difference between the top and the 
bottom of the layer i.e. A6 (t)AT st . Using such a temperature difference, it is then easy 
to define a time-dependent Marangoni Ma(t) = A6o(t)Ma or a time-dependent Rayleigh 
numbers Ra(t) — A9o(t)Ra. 

f Details are contained in Appendix D. 
J Details are contained in Appendix D. 
If In this case, we computed some additional results to cover a larger range of Biot numbers 
than the one given by Pearson in his article. 
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Let us compute the maximum of A9o(t) obtained during the time evolution for the 
Marangoni (rcsp. Rayleigh) case. This maximum is reached at time t c and leads to a new 
Marangoni number Ma(t c ) (resp. Rayleigh number Ra(t c )) which ca n be c ompared to the 
critica l values Masteady (resp. Rasteady) obtained bv lPearson ( 19581) (resp. lSparrow et aT 



(1963)) from the steady case. With the scalings used here, the critical Marangoni Ma c 
predicted from these steady results reads 

Masteady , D Ra-Steady 

Ma < * "a^T and Ra < * a^T (47) 

Such estimations have been plotted on figure 2] and [U They are found to be close to our 
results for Gthres — 1- The same comment applies for the critical wavenumbers except 
for the Benard-Marangoni case at high Biot numbers. For critical times, however, the 
steady-state approximatio n differs from o ur r esults. 

When using the results by IPearson ( 19581 ) or Sparrow et al. ( 1963^ . we are clearly using 



the normal mode results obtained for a linear temperature field in z on the whole thick- 
ness, which is a rough approximation especially at high Biot number. We can go even 
further and compare the non-normal mode results within the frozen-time approximation. 
In this latter approximation, a stability analysis in terms of normal modes is performed 
at the each time t. The "steady" base flow is assumed to be the temperature field 6q(z, t) 
computed at this specific time t. When the flow is stable within this frozen-time approx- 
imation for each time t, it will be assumed stable. When the control parameter reaches 
a critical value (here Ma c or Ra c ), there exists a unique critical time t c for which the 
frozen time state Oo(z,t c ) possesses a marginal eigenvector with a critical wavenumber. 
In the present problem we have determined the critical parameters for neutral conditions 
in the frozen-time case by the Lapack routine GGEV for generalized eigenvalue problems 
in combination with a Chebyshev collocation method. As can be seen in figures H] and 
results for the thresholds and the critical wavenumbers using the non-normal approach 
(Gthres = 1) an d the frozen-time approximation are close. However the prediction of crit- 
ical times t c still differs Q- The frozen-time approximation apparently provides a bound 
of order unity for t c because the available temperature difference A9$ (t) attains its max- 
imum as soon as the thermal boundary layer of the basic temperature distribution has 
reached the bottom of the layer. For small Bi, A0o(t) remains fairly constant for larger 
times, and the instability can develop on this quasi-steady background over fairly long 
times. This observation can explain the apparently unbounded growth of t c with Bi for 
Bi — ► in the non- normal analysis. 

4.5. Results for finite Prandtl numbers 

In this section, we focus on finite Prandtl numbers and more specifically on the role of 
initial perturbations on the transition zone estimation. In this case, the velocity field is 
not slaved to the temperature field so that perturbations in velocity can be considered as 
well as perturbations in temperature. We only discuss the curves for the pure Marangoni 
case (Ra = 0) (Figure [TJ but similar results apply to the pure Rayleigh case (Ma = 0). 



f In the frozen-time approximation, t c corresponds to a time when the quasi-static growth 
rate a(t) becomes zero. This critical time is of a different nature than the critical time for the 
non-normal approach. The latter characterizes the perturbation evolution over the time inter- 
val from t — up to t c . In zeroth order WKB theory, the critical time t c for the non- normal 
approach would be determined by f ' c a(t)dt — ln(Gth res ). Suitable modifications of the frozen- 
time analysis based on this observation should therefore lead to closer agreement regarding the 
critical times. We do not pursue this issue further in the present work. 
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Since the frozen time approximation seems to be valid for this unsteady problem and an 
exchange of stability in a normal mode is not affected by the Prandtl number, the effect 
of this latter number should not be significant. Indeed, for temperature perturbations, 
the critical Marangoni Ma c (Bi, Pr) is not affected when infinite Pr number case is 
compared to Pr — 10. For the optimal time t c (Bi, Pr) and wavenumber k c (Bi, Pr), 
the same conclusions apply. For velocity perturbations, the curves differ but not in a 
drastic way taking into account that perturbations are of a completely different nature 
compared with the infinite-Prandtl-number problem. Perturbing the temperature is more 
efficient than perturbing the velocity. Indeed critical Marangoni numbers are smaller for 
temperature perturbations, but the difference is of the same order as the one obtained by 
changing the threshold value from 1 to 100. Actually, the "blurriness" of this transient 
problem induced by the choice of the threshold values and perturbation types is not very 
broad and does not modify the order of magnitude of the critical numbers if one excepts 
the time t c . Finally, let us note that non -linear direct simulations have also been made 
to solve this problem (|Touazi et al. I (l2009h . showing very good agreement with the linear 
results presented here. 



5. Comparison with experiments 

The present section is devoted to the comparison between results ob tained from the 



optim al mode calculations and from an experimental work described in iToussaint et 



( 2008) , where transient Rayleigh-Benard-Marangoni convection is generated by drying a 
polymer solution of Poly Isobutylene- Toluene at ambient temperature. In the experiments, 
buoyancy and Marangoni effects are equally present. It is the evaporation of the solvent 
i.e. toluene which cools the upper surface by latent heat. During the experiments, the 
following parameters are kept constant: 



k = 0.97x10~ 7 m 2 S -\ A = 0.142 W K^ 1 m" 1 , p = 865 kgm' 3 , 
Q' = 1.07xl0~ 3 K~ x , cr = 28xl0~ 3 N m~\ 7 = 1.19 x 10~ 4 N K~ l m~~\ (5.1) 
L = 3.96 X 10 5 J kg' 1 , H th = 28 W K~ x m~ 2 , AT st = 4.8 K. 

Different thicknesses d and dynamic viscosities \x are considered, d is varied from 0.3 mm 
to 23.5 mm while dynamic viscosity \i is set to a value in the range [0.55 mPa s , 2100 mPa s] 
by monitoring the initial polymer concentration. These data allow us to estimate the 
crispation, the Galileo and the Peclet numbers for each experiment. We get the following 
bounds for these numbers: 

10~ 7 < CY < 10~ 3 , 2 x 10 2 Ga < 3 x 10 8 , 10~ 3 < Pe < 0.1 (5.2) 

As a consequence, the assumption of planar free surface and constant thickness layer are 
justified. The other relevant non-dimensional numbers vary in the following range : 

0.06 < Bi < 5; 6.6 < Pr < 2.5xl0 4 ; 20 < Ma 1.2x 10 5 ; 1.3 < Ra < 1.4x 10 6 (5.3) 

The comparison is displayed using the critical dynamic viscosity fi c (d) as a function of 
thickness d (figure [5]). In the plane (d,fi), each point corresponds to a given parameter 
set (Ra, Ma, Pr, Bi). The four critical curves correspond to two thresholds (Gthres = 1 
and Gthres = 100) and two types of initial pertubation (temperature and velocity). The 
theoretical critical curves divide the experimental points corresponding to regions of 
convection or pure conduction in a satisfactory manner. The temperature perturbation 
critical curve is above the velocity one. This analysis in the thickness/viscosity plane 



14 



F. Doumenc, T. Boeck, B. Guerrier and M. Rossi 



shows again that the bandwith of uncertainty due to the choice of threshold Gthres and 
perturbation types is not very broad and does not modify the order of magnitude of the 
critical thickness. 



6. Conclusion. 

This paper presents a novel linear stability analysis of an unsteady base state within 
the general conceptual framework of amplification theory. The non-normal approach is 
used, which possesses the advantage over more classical methods to solve the transient 
problem in a mathematically rigorous way. In turn, this allows one to test other ap- 
proximations. Here we have applied this approach for the first time to characterize the 
stability of a transient Rayleigh-Benard-Marangoni problem in an horizontal fluid layer 
suddenly cooled from above. It provides the upper limit of the energy amplification that 
a disturbance of wavenumber k can reach at time t. This quantity reaches a maximum at 
time t — t opt (Ma, Ra, Bi, Pr), for a specific optimal wavenumber k opt (Ma, Ra, Bi, Pr) 
and a specific perturbation structure in z. These latter two quantities play the role of the 
most amplified wavenumber and most amplified mode for the standard steady analysis. 
A "stability" diagram in the space (Bi, Pr, Ma, Ra) has been determined for the pure 
Marangoni and the pure Rayleigh problem by the non-normal approach. Note that the 
marginal conditions used to determine the stability curve was obtained by setting the 
optimal amplification equal to 1 or 100. The critical time and critical wawenumber were 
evaluated for this marginal conditions. Critical Marangoni and Rayleigh numbers exhibit 
a strong dependency on the Biot number and a weak sensitivity to Prandtl number vari- 
ations in the range Pr ^ 1 . Scaling exponents for critical Rayleigh or critical Marangoni 
versus Biot numbers have been found numerically and confirmed by scaling analysis in the 
limit of very small and very large Biot numbers. Comparison of the non- normal approach 
with the frozen-time approximation (a classical quasi-static approach) shows similar re- 
sults for the critical Marangoni or Rayleigh numbers and the critical wave numbers. 
Moreover, the "blurriness" inherent in any transient problem was analyzed as a function 
of the amplification threshold values and perturbation types. It has been shown that the 
transition region is thin compared to the large domain of Rayleigh or Marangoni numbers 
covered by the analysis. 

Finally, a comparison with experimental results has been performed, where convection is 
induced by solvent evaporation during the drying of polymer solution. A good agreement 
was indeed found between the present theoretical study and experimental observations. 
The method was developed in this paper for the cooling of a fluid induced by solvent 
evaporation but could easily be extended to other transient problems. 

TB and MR acknowledge financial support from the Deutsche Forschungsgemeinschaft 
in the framework of the Emmy-Noether Program (grant Bo 1668/2). 

Appendix A. Basic State 

We obtain here the scaling laws given in section 3 which provide the evolution of a purely 
conductive basic state. It is recalled that the basic temperature field is initialy uniform 
i.e. (z, t = 0) = 0. and satisfies a heat equation 



dt 



d 2 6 Q 



(Al) 



dz 2 
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with the boundary conditions 

d z 9 = at z=0, d z o + Bi0 o + Bi = O at z=l. (A 2) 

Due to the evaporation, a cooled layer of characteristic thickness <$(£) develops from the 
upper surface. Equation (|A 1[) provides the standard estimate 

5 (t) ~ min(Vt, 1). (A3) 

In the following, we determine the scaling laws for the two extreme cases Bi -C 1 and 
1 < Bi. 

,4 J Case Bi < 1 

Since 0o(z, i) is initially zero, it remains small during a period of time t < ti (the time 
Ti is determined below and shown to be much larger than 1). Diffusion in the liquid and 
evaporation terms thus dominate in the free-surface boundary condition (|A 2j) . Such a 
balance can be expressed in terms of order of magnitude as follows 

One thus obtains using equation (|A 3j) 

A6» - BiVi, 5 (t) ~ Vt for < Vt < 1, (A 5) 

During this time interval, the cooling layer has not reached the bottom at z = so that 
temperature field 9q reads in terms of orders of magnitude 

|0 O | ~ A0 O ~ BiVi for 0<Vi<l. (A 6) 

Thereafter the thickness remains constant 6o(i) ~ 1 and the following condition holds 

A0 o (t) ~ Bi, S {t) - 1 for 1 < t < ri. (A 7) 

During this latter time interval, the heat equation (|A ip can be used with scaling (|A 7p 
to get the temperature field 0q in terms of orders of magnitude 

\0 o \~Bit, for l<*<n. (A 8) 

These equations are valid if evaporation dominates heat transfer in the gas phase. This 
requires that \8q\ <C 1 and determines the value t± ~ Bi^ 1 . For t± < t, the temperature 
field #o relaxes towards the steady uniform temperature equal to 0q(z, t) = — 1. Since the 
temperature gradient is equal to zero in that regime, all the energy due to evaporation 
is transfered by convection in the gas phase. 

B) Case 

For small times, the condition |6*o| <C 1 holds and an analysis similar to the one performed 
for the case Bi <C 1 is valid leading to 

A6> ~ BiVi, 6 (t)~Vi for < t < r 2 , where T 2 ~Bf 2 . (A 9) 

The value of Ti is obtained by determining the time when A#o ~ 1. At that time, the 
heat flux in the gas phase becomes of the same order of the evaporation. During this 
time interval, the cooling layer has not reached the bottom at z — so that temperature 
field £>o reads in terms of orders of magnitude 

1 0o | - A6» - BiVi for < t < Br 2 . (A 10) 
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For Bi 2 < t, a new regime begins where the surface temperature remains constant 
\6q(z — 1, t) \ ~ 1 while the cooled layer thickness keeps increasing 

A6> ~ 1, 5 ~ V7 for time Br 2 < t < 1. (All) 

This regime ends when Jo ~ 1 at time t ~ 1. Thereafter the temperature decreases in 
the whole layer thickness to reach the steady state regime 9q(z, t) = — 1. 
The case Bi ~ 1 is the limiting case of the two previous ones. The equation (|A5p is valid 
until t ~ 1, when the cooled layer thickness and the surface temperature both reach their 
extremum. Thereafter the temperature decreases in the whole layer thickness, to reach 
the steady state. 



Appendix B. Obtaining the Adjoint Equations 

Let us denote by qj(z,t),j = 1, ... ,4, the components of the vector field 

{u(z,t),w(z,t),6(z,t),p(z,t)). 

To find the maximum amplification at a given time t\, we maximize the perturbation 
norm E{q{t{)) 

3 

E(q(ti)) = ^2 Q J «&(*,ti)g+(2,ti)cU! (Bl) 

3=1 

at time ii with respect to the set of all possible initial perturbations q(0) such that 
^(9(Q)) — 1- It is recalled that the integration is performed over the entire layer depth 
and the superscript + denotes complex conjugation. Coefficients Cj are weight coefficients 
chosen to put emphasis on temperature or velocity acoording to the case considered. To 
analyse the initial perturbation in velocity, the kinetic energy norm Ey is used and one 
takes G\ = G% = 1 and C3 = 0. To analyse the initial perturbation in temperature, the 
temperature norm Et is used and C\ — C2 = and C3 = 1. 

The variation SE(q(ti)) with respect to a variation Sq(Q) of the initial perturbation is to 
be evaluated. This computation cannot be performed in a straightforward manner since 
the energy E(q(t\)) can be explicitly written in terms of q(ti) but only implicitely in 
terms of q(0). It is known via several constraints: normalization of q(0) and time integra- 
tion over the interval [0, £1] of equations (3.8)-(3.10). These dynamic equations relating 
q(0) to q(ti), are formally written here as Fj (q) — 0, j = 1 . . .4. This optimization with 
constraints necessitates the introduction of Lagrangian multipliers, the so-called adjoint 
fields q(t) = (u(z,t),w(z,t),d(z,t),p(z,t)). 

More specifically, a Lagrangian function L is defined, which depends on direct q(t) and 
adjoint q(t) variables over the interval [0, ti], and a normalization scalar Sq: 

4 -ti 

L( 9) g, So ,*i) =E(q(t 1 )) - s o (E(q(0)) - 1) - ^ / alt ({Fj (q(t j) , qj (t)) + c.c.) (B2) 



where c.c. means complex conjugate and (•, •) stands for the scalar product 

(01,02)= J ai(z)&2(z)dz. (B 3) 

When q(t) satisfies the constraints (direct problem plus normalization at t = 0), all terms 
but the first one on the r.h.s. of equation (|B 2j) are zero and, by consequence, L = E and 
SL = SE. At this stage, the adjoint variables and the quantity sq are left unspecified. 
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Formally the variation SL reads as 

SL = Y,Cj([ qt(z,t 1 )Sq j (z,t 1 )dz-s J q+(z, 0)Sq, ;(z, 0)dz\ (B4) 
j=i ^ / 

- J2 / dt[(5F 3 (q(t))A 3 (t)) + {F S (qm <%(*)>] + c - c - 

The expression (Fj(q(t)),Sqj(t)} in equation (|B 4|) is zero if the governing equations 
Fj(o) = are satisfied during the time interval [0, ti]. The main idea then amounts to 
rewriting quantity (5Fj(q(t)), qj(t)) in terms of Sqk(t). This is done by integrating by 
parts in space or time. After some tedious algebra, the following identity 

4 r h 4 tl 



]T / 1 dt(SF s {q(t)),gj(jt)) - Et / 1 dt(^(«(t),«),*g(t)>] + (B5) 

2 



3 = 1 

+ J (gfatJSqafatJdz- f qt(z,0)6q 3 (z,0)dz + B(Sq,q) 

can be established, where Fj is an expression containing spatial or time derivatives of q. 
Note that the second, third and fourth r.h.s terms originate from integrations by parts 
of time derivatives in equations (3.8)-(3.10) and terms B(Sq,q) are generated from the 
boundary terms resulting from integrations by parts of spatial derivatives. These latter 
terms hence involve only quantities Sq and q evaluated at the boundaries z = and 

3 = 1. 

At this stage, the freedom of the Lagrangian multipliers can be used to impose some 
added constraints on the adjoints fields q, namely: (i) equations Fj(q(t)) = 0,j = 1...4, 
which are similar to Fj(q(i)) for q and define the evolution equations (3.15)-(3.17); and 
(ii) boundary conditions B(5q, q) = 0, which are the counterpart of boundary conditions 
(3.10)-(3.12) on q and define the boundary conditions (3.18)-(3.19) for q. This new 
system can be simulated as the direct problem. It is easily seen that the adjoint system 
(3.15)-(3.19) must be integrated backwards in time. When it is satisfied, the variation 
SL reads 

SL = J2[J (C jq f(zM) - -^£+(2, *!))<%(*, ti)d3 (B6) 



(s C jq +(z, 0) - 0))<%(3, 0)dz 



(s C 3 q^(z,Q) - q+(z,0))S qj (z,0)dz + c.c. 



Two relations can be still imposed, a first one at time t = t\ which relates qj{z,t{) and 
qj(z, ti) and a second one at time t = which relates qj(z, 0) and </j(z,0). These two 
constraints are satisfied so that SL = and are defined precisely below according to the 
norm and Prandtl number. The condition SL — means that an optimal perturbation is 
attained. However this process should be self-consistent : one uses the iteration procedure 
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(3.20). When the iterative process has converged, an initial optimal perturbation for time 
ti is found. 

A) Finite Prandtl and zero initial temperature perturbations 

If one considers only initial perturbations in velocity field so that variation of temperature 
field Sq 3 (z, 0) is zero, it is consistent to use the norm Ey, i.e., C± = Gi = 1 and C 3 = 0. 
Equation (|B 6p then naturally leads to the relation 



Qj{z,h) = Prqj(z,ti), j = l,2 
at time t = t\ and the relation 



qj(z,0) = 



1 



Prs 



9,(2,0), 3 = 1,2 



$3(2, h) = 



93(^,0) =0 



(B7) 



(BJ 



at time t = 0, where so is chosen such that the normalization condition E(q(0)) = 1 is 
satisfied. 

B) Finite Prandtl and zero initial velocity perturbations 

When considering only initial perturbations in temperature field so that 5q%(z, 0) = 
Sq2(z,0) = 0, it is consistent to use the norm Et, i.e., C\ = C2 = and C3 = 1. 
Equation (|B 6p then leads to the relation 



q j (z,t 1 ) = j = l,2; q 3 (z, h) = q 3 (z, h) 



at time t = t\ and 



gj (z,0)=0, i = 1, 2 



93(^,0) = —q 3 (z,0) 
so 



(B9) 
(B10) 



at time t = so that the normalization is satisfied. 
C) Infinite Prandtl number 

For the infinite Prandtl number, the norm Et is chosen since the velocity is slaved to 
the temperature field in that instance, hence C\ = C2 = and C3 = 1. The equations 
are then identical to the previous case except that only the equation for temperature 
appears, i.e., 



and at time t = 



<&0Mi) = Q3(z,ti), 



1 



93(2,0) = —93(21,0) 
so 

so that the normalization is satisfied. 



(Bll) 
(B12) 



Appendix C. Numerical Method. 

For the numerical solution, the direct and adjoint equations are reformulated as fourth- 
order problem in a streamfunction-like approach. The incompressibility constraint is 
thereby satisfied automatically, and the pressure and horizontal velocity are eliminated 
from the equations. For finite Prandtl number, the direct equations take the form 



J_d_ 



V 



d6o 

dz 



dz 2 
dz 2 
dz 2 



fj - Rak 2 6, 



w, 



(Cl) 



(C2) 



(C3) 
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The boundary conditions are 



10 



dw 
~dz~ 



d9_ 

dz 



— = — = at z = 



and 



w = 0, fj + k 2 MaO 
The adjoint fields satisfy the system 

1 d . 



Pr dr 



n = 



v 



= o, 

dz 2 

cP_ 

dz 2 
dz 2 



89 
dz 

-k 2 
-k 2 



BiO = at z = 1. 



220*0 



fj-k z e 



+ Raw 



with the boundary conditions 



— = — = atz = 

oz oz 



(C4) 

(C5) 

(C6) 
(C7) 
(C8) 

(C9) 



and 



w 



fj = 0, 



BiO + Ma 



dw 



at z = 1. 



(C10) 
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These equations are discretized in time with a backward Euler method for the diffusive 
terms. The product term with the basic temperature profile is treated with the explicit 
Euler method. For the direct problem the solution at the new time level n + 1 is obtained 
by solving the following equations in sequence: 

r d 2 



dz 2 



k 



dz 2 



1 

" At 
1 



PrAt 



nn+l 



p.n+1 



At 

fj" 



PrAt 



dz ' 

Rak 2 6 n+1 , 



dz 2 



- k< 



(Cll) 
(C12) 
(C13) 



The boundary equations for fj are given in terms of w. To satisfy them, the solution of 
the second and third equation is represented by the linear combination 



V 



VP - 

Wp 



Xfjl + [lf]2, 
- \Wl + /J,U>2, 



(C14) 
(C15) 



where the solution with subscript P is a particular solution of the r)-equation (jC 12|) with 
fjp = on the boundaries z — and z = 1. The functions with subscripts 1 and 2 satisfy 
the homogeneous 77-equation with zero right hand side and two linearly independent 
boundary conditions, which we choose as 



f)x{z = 1) = fji(z = 0) = 1, 
fj 2 (z = 1) = -fj 2 (z = 0) = 1. 



(C16) 
(C17) 



The boundary conditions dw/dz = at z = and fj + k 2 Ma9 = at z = 1 determine 
the coefficients A and fi in the linear combination. We note that the functions wp, w\ 
and W2 satisfy zero boundary conditions at z = and z = 1. 
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The adjoint solution at the new time level n + 1 is likewise found by solving the equations 

d 2 t3 _ _J_1 ~n + X = _J[_ _ tffrW 

PrAr ' PrAr dz ' 



dz 2 



dz 2 



w n+1 



r) 2 1 

2. i- 2 — 

Dz 2 At 



qn+l _ 



Raw 



~n+l 



(C18) 

(C19) 
(C20) 



in the given sequence. The solution for w n+1 and fj n+1 must again be represented as 
a linear combination with auxiliary functions satisfying the homogeneous ry-equation in 
order to satisfy the boundary conditions. 

Discretization in sp ace is realized with an expansion in Chebyshev polynomials (see 
Canuto et all (1988)). Product terms with the perturbation and the basic state are cal- 
culated in physical space by a fast cosine transform. The Helmholtz equations for the 
variables 77, w, 9 and the adjoint variables fj, w, 9 reduce to essentially tridiagonal linear 
systems. The boundary conditions are treated with the tau method, which produces two 
filled rows in the matrix representation. The limit of infinite Prandtl number requires no 
changes in the solution procedure. 

The basic temperature profile is computed with the same numerical method as the per- 
turbations, i.e. using the backward-Euler method and a Chebyshev polynomial expansion 
with the same time step and number of polynomials. The entire field 9q(z, t) is stored in 
an array in order to speed up the backward integration of the adjoint equations. 
The code was tested for infinite Prandtl number with a stationary basic temperature 
profile. It was verified that exponential growth of the optimal perturbations appeared at 
the proper thr eshold values of Ma ss 79.6 for pure Marangoni convection with Bi = 
(jPearsonlll958T ) and for Ra sa 1100 for pu re Rayleigh convection with fixed temperature 
on the free surface (| Chandr asekharl 1 1 96 if ) For this verification, the boundary condition 
at the bottom was changed to constant temperature. 



Appendix D. Analysis of Stability for the Marangoni case 

D.l. The approach for the Marangoni case (Ra = 0). 

This section presents an approach valid for infinite Prandtl number, which evaluates the 
evolution in terms of orders of magnitude. It is based on two hypotheses which make the 
analysis tractable. First the initial perturbation of wavenumber k in the x direction is 
only a temperature perturbation i.e. u(z,t = 0) = w(z,t = 0) = and the temperature 
perturbation 9{z, t = 0) is uniform along the z-direction. Second, the flow is supposed 
unstable i.e. convection sets in, if there exists a time and a region in the flow in which the 
advection term in equation (3.10) becomes greater or equal to the two diffusion terms 
which tend to damp the initial perturbation. 

Practically, one first determines the orders of magnitude for temperature perturbations 
when the system is in a stable regime or near the critical curve, i.e., when the term corre- 
sponding to advected heat transfer w^f in equation (3.10) can be neglected, according 
to the second hypothesis. Thereafter one computes the order of magnitude of the term 
k 2 9 i.e., the diffusion in the cc-direction, and of the term 5-4, i.e., the diffusion in the 
z-direction, corresponding to the stable regime. This is done in subsection ID. 21 On the 
other hand, the order of magnitude for velocity w is found in subsection ID. 31 as well as 
the corresponding advection term, w^-. 
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Since the velocity is computed using the temperature perturbation field estimated for 
the stable configuration, this approach is consistent only if the advection term remains 
much smaller than one of the diffusion terms. The sets of parameters Ma, Bi, Pr that 
give consistent results for each time and mode are considered in the "stable" domain. 
Otherwise, if there exists a time and a mode of wavenumber k such that the advection 
term is larger in order of magnitude than the two diffusion terms, the corresponding set 
of parameters is associated with a situation where convection sets in (subsection ID.4[) . 
The scaling laws for critical parameters are then derived by solving the resulting set of 
inequalities (subsection ID. 5|) . 

D.2. Scaling Analysis for the Temperature Perturbation Field 

As mentioned in the previous paragraph, we determine the orders of magnitude for 
temperature perturbations by neglecting the advected heat transfer w^- in equation 
(3.10). One thus obtains 

00 
~di 

The temperature perturbation field also satisfies the boundary conditions 



d 2 
dz 2 



= (D 1) 



09 

— +Bi6 = 0,&tz = l, (D2) 
oz 

f)0 

— = at z = 0. (D 3) 

oz 

The cooling due to evaporation imposes that a thermal boundary layer for the tempera- 
ture perturbation 6(z, t) develops. One can easily check that the solution 

6(z, t) = o (l + (z, t)) exp{-k 2 t) (D 4) 

satisfies the above system and the condition of uniformity at t = 0. Note that ao is 
simply the initial amplitude of the temperature perturbation which is taken to be equal 
to one in the sequel. From the above equation, it is readily seen that the thickness of the 
thermal boundary layer for the perturbation field 9{t) is equal to 5o(t) and that the scale 
of variation in the ^-direction of the perturbed field 0{z, t) denoted by A0 satisfies 

A9(t) ~ A0 o {t) exp(-fc 2 t) (D 5) 

The scale S of the perturbed field 0(z = l,t) on the surface satisfies according to 
equation (|D 2[) 

§ 1 A§ ( f ) 
s Bi 6 {t) 

Using the results of Annex A, it is straightforward to find the following estimates 
For Bi < 1 : 



(D6) 



A9(t) - BiVtexp(~k 2 t), 6 S - exp(~k 2 t), - exp(-k 2 t) for 0<Vi<l, (D 7) 

A0(t) - Biexp(-k 2 t), S ~exp(-fc 2 i), 6 - exp(-fc 2 t) for 1 < t < Bi' 1 . (D 8) 
For l<Bi: 

A0(t) - BiVtexp(-k 2 t), S - exp(-fc 2 i), 6 - exp(-fc 2 <) for < t < Br 2 , (D 9) 
A9(t) ~ exp(-fc 2 <), S ~ CX P(~^) | §^ cxp(-k 2 t), for Bi~ 2 < t < 1. (D 10) 

Biy't 
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D.3. Scalings for the Velocity Field in the Benard-Marangoni problem (Ra = 0). 

The equation of the vorticity field can be easily deduced from equations (3.8) and (3.9). 
Denoting by uj the y-componcnt of vorticity, we obtain the diffusion equation 

—d t uj = d 2 z Cj -k 2 6j (D 11) 

Pr 

For infinite Prandtl number, this equations simplifies 

d 2 z uj -k 2 u} = (D 12) 

The vorticity is slaved to the temperature evolution via the boundary condition at the 
free surface given by equation (3.11) : 

Cu + ikMa6 = at z = 1. (D 13) 

Equation (|D 12[) plus the forcing (|D 1 3[) defines an hydrodynamic boundary layer Sh- It 
is easily seen that the proper scaling reads 

S H ~min(l/fc, 1) (D14) 

The hydrodynamic boundary layer either reaches the bottom, i.e., Sh ~ 1, or the diffusion 
term along x becomes of the same order of the diffusion term along z and 5^ ~ 
From equation (3.11), a relation between 6 S and the order of magnitude of velocity u can 
be found : 

u~k6 H Ma6 s {t). (D15) 
The order of magnitude of the vertical component w of velocity is provided via mass 
conservation 



w ~ kS H u = (kS H ) 2 Ma 6 s (t). (D 16) 

D.4. Condition for the onset of convection for the Marangoni flow (Ra — 0) 

To describe the time evolution of perturbations, one must distinguish two regions along 
the z direction i.e. inside and outside the thermal boundary layer. Outside the layer 
(^o < 1~ z l)i the advection term in equation (3.10) is zero since the basic temperature 
field 6$(z,t) vanishes : hence diffusion dominates and perturbations are always damped. 
Instability thus only arises inside the thermal layer (0 < 1 — z < <5o)- 
To determine the onset of instability, one first compares the order of magnitude of the 
advection and the diffusion along the x-direction in the thermal layer 

^f 90 , kH (D17) 
do 

where Wth denotes the order of magnitude of the vertical velocity in the thermal boundary 
layer, and second the order of magnitude of the advection term and of the diffusion term 
in the z-direction 

™thM AO 

^o~^ (D18) 
The existence of a convection onset thus implies that there exists a time and a mode of 
wavenumber k for which the two conditions 

m h (^)>kH and ^(^)>^ (D19) 
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hold. We need quantity wth since it explicitely appears in the above inequalities. Two 
possibilities should be considered at each time : So(t) < Sh or Sh < So(t). In the first 
case, the thermal layer is included in the hydrodynamic layer and one may use the scaling 
w t h ~ j 2- ^ ln the second case wth ~ w. Using equation (|D 1 6[) . this implies that the 
scaling Wth is such that 

w th ~ min (l, w = min (l, |0 (kd H ) 2 Ma 9 S , (D 20) 

It is now possible to rewrite inequalities (|D 19|) as 

min f i ; *±) 5 2 H Ma 6 S (^) > 9 and min ( 1, ^] (kS H ) 2 Ma 9 S A9 > ^ (D 21) 
\ S H J o \ o H J d 

D.5. Derivation of scaling laws 

One must now introduce the various expressions previously obtained for A9q, So, 9, A9, 
6 S , Sh inside instability conditions ID 211 The expressions for So and A9q are obtained in 
section (jXJ, the expressions for 9, A9 and 9 S in section (|D.2|) . the expression for Sh in 
section (|D.3|) . To ease the discussion, three separate cases are studied : 
A) Bi < 1 , B) 1< Bi and t < Bi~ 2 , C) 1< Bi and Br 2 < t. 

A) Case Bi < 1 

First let us recall from Annex A, that the following relations hold : 

5 (t) ~min(Vt,l) and ~ Bi. (D22) 

From section fD. 21 one easily verifies that the temperature perturbation field is such that 

A9 ~ BiS 9 s and 9 ~ 9 S (D 23) 

Using equations (|D 22|1 and (|D 23|) . condition (|D 21[) can be transformed into 

min f 1, p-^j Ma Bi S 2 H > 1 and min f 1, |^ 5 Ma fc 2 5| > 1 (D 24) 

Note that only the period t < l/Si should be considered here since, for 1/Bi < t, the 
basic state has relaxed to a uniform temperature. 

To ease the discussion, two separate cases must be considered for the wavenumber k, 

namely k < 1 and 1 < k. 

•k<l 

In that instance, Sh ~ min(l/fc, 1) = 1 (see equation (|D 14]) ) leading to the equality 
min(l, |i) ~ 5 . Condition (IET24)) reads 

S Ma Bi > 1 and S 2 Ma k 2 > 1 (D 25) 

The smallest Marangoni number i.e. the critical Marangoni number which satisfies such 
inequalities, is obtained for So(t) ~ 1 i.e. for 1 < t. Condition (|D 25|) becomes 

Ma > 4t and fc 2 > (D 26) 

~ Bi ~ Ma v ' 

From the above onditions, one easily gets the critical value 

Ma c ~l/Bi, with ^Bi < k c < 1 and l<t c <l/Bi. (D 27) 



24 F. Doumenc. T. Boeck, B. Guerrier and M. Rossi 

In that instance, Sh ~ min(l/fc,l) = l/k. Since the two functions min(l,fc(5o) and 
min(l, kSo)So are both increasing functions of So when So < 1/fc, the critical Marangoni 
must be obtained when l/k < So for which min(l, kSo) = 1. Conditions (|D 24[) become 

Ma Bi >k 2 and Ma S [t) >1 (D 28) 

A straigthforward discussion directly leads to the conditions 

Ma c ~l/Bi and fc c ~ 1 and \<t c <l/Bi. (D 29) 

which is a limiting case of condition (|D 27[) . As a consequence, the critical conditions in 
the case Bi <g; 1 corresponds to condition (|D 27[) . 

B) Case 1 « Bi and t< Bir 2 . 

From results obtained on the basic flow, it is easily seen that relations (|D 22[) . (|D 23|) 
and thus (|D 24p still hold. Moreover note that the largest value of So(t) is obtained at 
the largest time t ~ Bi~ 2 : So(Bi~ 2 ) ~ Bi^ 1 . One must consider the three cases fc < 1, 
^^k<Bi and < fc. 

• A; < 1 

In that case, 5# ~ min(l/fc,l) = 1 and condition (|D 25|1 is again satisfied. The critical 
Marangoni number, is obtained for the largest possible value of 8o(t) i.e. 5o(Bi~ 2 ) = 
Bi -1 . Condition (|D 25|) now reads 

Ma > 1 and Ma k 2 > Bi 2 (D 30) 

A straigthforward discussion directly leads to the critical conditions for instability 

Ma c ~Bi 2 with k c ~ 1 and t c ~ Br 2 (D31) 

• 1 < k < Bi 

In that instance, Sh ~ min(l/fc, 1) = l/k. Conditions (|D 24[) become 

min(l, jfe5 (*)) fc -2 Ma Bi > 1 and min(l, kS (t)) S (t) Ma > 1 (D32) 

The critical Marangoni number, is obtained for the largest possible value of So(t) obtained 
at the largest time i.e. t = Bi~ 2 for which So(Bi~ 2 ) = Bi -1 . As a consequence 

min(l, k/Bi) Ar 2 Ma Bi > 1 and mm(l,k/Bi) Br 1 Ma > 1 (D 33) 

For this wavenumber interval, min(l, k/Bi) = k/Bi and equation (|D 33[) becomes 

k- 1 Ma > 1 and k Br 2 Ma > 1 (D34) 

This leads to the critical conditions for instability 

Ma c ~ Bi with k c ~ _Bi and t c ~ i?i -2 (D 35) 

• Bi < k 

In that case, Sh ~ min(l/fc, 1) = 1/fc and conditions (|D 32[) are verified. Again the critical 
Marangoni number, is obtained for the largest time i.e. t = Bi~ 2 corresponding to the 
largest possible value of So(t). Moreover, for this wavenumber interval, min(l, kSo(t)) = 1 
and equation (|D 32j) becomes 

k' 2 Ma Bi >1 and Bi~ l Ma >1 (D36) 

This leads to the same critical conditions (|D 35|) for instability: 



Ma c ~ i?i with fc c ~ _Bi and i c ~ Bi 



(D 37) 
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C) Case and Bi~ 2 <t<l 

First let us recall from results on the basic flow that the following relation holds: 

5 (t)~^/i and A0 O (*)~1 (D 38) 

From section uT2l one easily verifies that the temperature perturbation field is such that 

A6 ~ BiS e s and ~ A9 (D 39) 

Using equations (|D 38|) . (|D 39|) . conditions (|D 21[) can be transformed into 

min fl, ^ 2 £| Ma Br 1 > 1 and min (l, 5% k 2 Ma Br 1 > 1 (D40) 

At this stage, two possibilites should be considered : k < 1 and 1 < k. 
•k<l 

In that case, Sh ~ min(l/fc, 1) = 1 and min(l, j 2 -) ~ So- Conditions (|D 40|) become 

Sq 1 MaBf^l and <5o fc 2 Ma Bi^ 1 > 1. (D41) 

By multiplying both conditions, one gets: 

Ma > Bik- 1 . (D42) 

This implies that k c ~ 1 and Mo c ~ Introducing the latter two equalities back into 
equation (|D 41j) provides 5q = 1 i.e. t c ~ 1. Finally the critical conditions can be written 

as 

Ma c ~Bi with fc c ~ 1 and t c ~ 1. (D43) 

• 1 ^ 

In that instance, 5// ~ min(l/fc,l) = 1/fc. If one introduces the new variable £ = kSo, 
conditions (|D 40|l read 



Ma > BiF(£) and Ma > BiG(£) (D 44) 

in which 

F(0 ee ^ 2 min(l,C) _1 with G(0=min(l,0 _1 - (D45) 
A straightforward analysis of these two functions shows that the critical Marangoni is 
reached for £ = 1 hence Ma c ~ Bi. Moreover, since 1/Bi < So < 1 in this time interval, 
a large bandwith of modes k are equivalent leading to the following critical conditions 
for instability: 

Ma c ~Bi with 1 < k c < Bi and t c ~ k~ 2 (D46) 

Finally, by taking the lowest Marangoni numbers of the conditions (|D 31|) - (|D 35|) - (|D 43jl - 
(|D 46[) , one deduces the true critical conditions for 1 <C Bi, namely 

Ma c ~Bi with 1 < k c < Bi and t c ~k~ 2 . (D47) 
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Figure 1. Basic temperature profile for different Biot numbers, (a-c) Temperature profile at 
times t — 1CP 2 , t — 1CP 1 and t = 1. (d) Temperature difference A6o(t) as a function of time t. 




Figure 2. The maximum energy amplification G(t;k; Ma, Ra, Bi, Pr) as a function of 
wavenumber k for three different times t — 0.2, t = 0.43 and t = 1. Parameters Ma = 300, 
i?a = 0, = 1, Pr = oo. 
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Figure 3. Isolines of the maximum amplification G ma x(Ma, Ra, Bi, Pr) in the plane 
{Bi, Ma) Parameters Ra — 0, Pr = oo. The values of the isolines are written on the figure. 





Figure 4. Infinite Prandtl number case and Ra = 0. Results are shown for two thresholds 
Gthres = 1 and Gthres = 100. Frozen-time and steady state results are presented for comparison 
(see text for details), (a) Critical Marangoni Ma c (Bi), (b) Critical wavenumber k c (Bi), (c) 
Critical time t c {Bi). 
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Figure 5. The optimal temperature perturbation 9(z,t) at time t = (dashed) and t = t c 
(solid) for Ra = 0, Pr = oo, G thres = 1 : (a) Bi = 0.01, Ma = Ma c = 8685 and k = k c = 0.74, 
(b) _Bi = 100, Ma = Ma c = 1658 and k = k c = 4.46. 
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Figure 6. Infinite Prandtl number case and Ma = 0. Results are shown for two thresholds 
Gthres = 1 and Gthres = 100. Frozen-time and steady state results are presented for comparison 
(see text for details), (a) Critical Rayleigh Ra c (Bi), (b) Critical wavenumber k c (Bi), (c) Critical 
time tc(Bi). 
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Figure 7. Infinite Prandtl ( Pr — oo) or finite Prandtl (Pr = 10) number cases. Temperature or 
velocity initial perturbation results with threshold Gthres = 1- (a) Critical Marangoni Ma c (Bi), 
(b) Critical wavenumber k c (Bi), (c) Critical time t c (Bi). 




initial thickness / 



Figure 8. Comparison of theoretical and experimental results in the plane (d, fj.) where d stands 
for the layer thickness and fj, the dynamic viscosity. The various curves displayed in solid lines 
correspond to two different thresholds Gthres = 1, Gthres = 100 and two different perturbation 
types (velocity and temperature). Experimental data are displayed by symbols. 
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